Aim: Detect pathways/genes that are differentially expressed in response to disease.

1 Background

1.1 Patient metadata

The 5 patient samples of interest are clinically diagnosed with mitochondrial disease and genetically diagnosed with mutations in the gene ATG7. Details below:

load(here::here("results/aberrant_expression/OUTRIDER/50_controls/gene_count_ods.rda"))

gene_count_ods <- gene_count_ods[, colData(gene_count_ods)[["samp_id_tidy"]] != "S2557"]

which_ATG7 <- which(colData(gene_count_ods)[["gene_name"]] == "ATG7")
which_ATG7 <- colData(gene_count_ods)[["samp_id_tidy"]][which_ATG7] %>% 
  as.character()

ATG7_patient_metadata <- colData(gene_count_ods)[which_ATG7,] %>% 
  as.tibble() %>% 
  dplyr::select(samp_id_tidy, gene_name, mutation_type, mutation_hgvs)

ATG7_patient_metadata$mutation_type <- ATG7_patient_metadata$mutation_type %>% lapply(FUN = str_c, collapse = ";") %>% 
  unlist()

ATG7_patient_metadata$mutation_hgvs <- ATG7_patient_metadata$mutation_hgvs %>% lapply(FUN = str_c, collapse = ";") %>% unlist()

ATG7_patient_metadata <- ATG7_patient_metadata %>% 
  mutate(family = c(2, 2, 3, 4, 1)) %>% 
  arrange(family) 

ATG7_patient_metadata %>% 
  data.table() 
ATG7_patient_metadata %>%
  write_delim(here::here("results/tables_figures/ATG7_patient_metadata.csv"), delim = ",")

1.3 OUTRIDER

  • OUTRIDER is a bioinformatics software that aims at detecting aberrantly expressed genes from RNA-seq data. It is specifically designed for a rare disease context (1 vs all) by correcting for co-variations that are not known a-priori using an autoencoder to to help distinguish outliers more readily. OUTRIDER then applies a negative binomial test to find outlier expressed genes.

1.4 DESeq2

  • DESeq2 aims at finding differentially expressed genes by modelling reads using a negative binomial distribution (similar to OUTRIDER), then using a binomial test to ascertain significance.

2 Method

Patient and control samples

RNA-sequencing was performed on a total of 58 individuals. 5 of these were ATG7 patients detailed in table X. In particular, samples were obtained from patients M1856.17, M0920.18, M0921.18, M1111.19 and M1716.19 from families 1, 2, 2, 3 and 4 respectively. The remaining 53 samples were used as the control samples consisting of a cohort of 46 patients diagnosed with mitochondrial disease, 7 patients with congenital muscular dystrophyand 2 unaffected individuals.

Culturing fibroblasts

Fibroblast cell lines were obtained from 58 individuals were obtained through collaboration with Dr. Haiyan Zhou and the Lily Consortium consisting of Dr. Charu Deshpande, Dr. Ines Barbosa, Prof. Joanna Poulton, Prof. Michael Simpson, Prof. Robert McFarland, Dr. Robert Pitceathly and Prof. Robert Taylor. DMEM supplemented with 10% Fetal Bovine Serum and 0.05 g/ml Uridine was used for culturing cells. Harvesting cells was performed by detaching cells using TrypLE Enzyme, washing cells once with washed with DPBS before storage at -80°C.

RNA-sequencing

RNA from fibroblast pellets was either extracted and sequenced at Eurofins Genomics or extracted through Bioxpedia, then sequenced at UCL Genomics. RNA integrity numbers (RIN) were measured using Agilent Technologies 2100 Bioanalyzer or Agilent 4200 Tapestation by Eurofins/Bioxpedia respectively. All RIN values exceeded 8.0. RNA library preparation was performed using the INVIEW Transcriptome Discover protocol from Eurofins. Specifically, the cDNA library was prepared using a random primed, strand specific, poly-A selection protocol. Paired-end sequencing was performed on illumina NovaSeq 6000 Sequencing System machines to an coverage of ~100 million pair-end reads per sample. Reads of 150-bp length were used to increase the proportion of reads with a gapped alignment, which represent splicing events.

RNA-sequencing alignment and quality control of patient samples

Pre-alignment quality control including adapter trimming, read filtering and base correction was performed using fastp, an all-in-one FASTQ preprocessor (v0.20.0) (CITATION - https://academic.oup.com/bioinformatics/article/34/17/i884/5093234). Reads were aligned using STAR 2-pass (v2.7.0) to the hg38 build of the reference human genome (hg38) using gene annotation from Ensembl v97 (CITATION - https://pubmed.ncbi.nlm.nih.gov/23104886/). 1st pass alignment was used to discover novel junctions, which were used as input to the 2nd pass to improve the sensitivity of junction detection. Via the option --outFilterMultimapNmax 1 reads were required to be uniquely mapped to only a single position in the genome. The minimum required overhang length of an annotated and unannotated junction was required to be 8 and 3 base pairs respectively. The output BAM files underwent post-alignment QC using RSeQC (CITATION - https://academic.oup.com/bioinformatics/article/28/11/1530/267467) with all samples passing quality control after manual assessment.

Detecting abberantly expressed genes and pathways

Gene count matrices were generated using RNA-SeQC (v2.3.4) with Ensembl v97 reference annotation matching the protocol used in https://github.com/broadinstitute/gtex-pipeline (CITATION - https://pubmed.ncbi.nlm.nih.gov/22539670/). OUTRIDER was used to find outlier expression of genes in each ATG7 patient, hence an all vs 1 experiment design was used, whereby each ATG7 patient were compared to the remaining 53 controls (CITATION - https://www.sciencedirect.com/science/article/pii/S0002929718304014?via%3Dihub). OUTRIDER was run with default setting and the autoencoder was fit using the optimal encoding dimension found to the be XX. P-values were corrected for multiple testing using the Benjamin Hochberg method and 0.05 was used as the signifance threshold. DESeq2 was applied to detect common genes or pathways that were dysregulated across all ATG7 patients (CITATION - https://genomebiology.biomedcentral.com/articles/10.1186/s13059-014-0550-8). Thus, all 5 ATG7 were grouped together as the case cohort and compared to the 53 control samples. The Benjamin Hochberg method was used for multiple test correction and genes with a p-value less than 0.05 and an absolute log-fold change greater than 1.5 were considered differentially expressed. gprofiler2 was then run to with differentially expressed genes as input, ranked by p-value to find dysregulated pathways (CITATION - https://academic.oup.com/nar/article/47/W1/W191/5486750).

Detecting abberantly expressed splicing events

dasper (https://github.com/dzhang32/dasper) was used to detect aberrant splicing events in each of the 5 ATG7 patients. Foremost, BAM files were converted into the BigWig format using megadepth (CITATION - https://github.com/ChristopherWilks/megadepth). Junctions outputted by STAR and BigWig files were inputted into the dasper pipeline. Foremost, dasper loads, filters junction data. Here, we required junctions to have at least 5 counts in at least 1 sample, a length between 20-1,000,000 base pairs and to not overlap any ENCODE blacklist regions (CITATION - https://www.nature.com/articles/s41598-019-45839-z). Junctions were then annotated using Ensembl v97 reference annotation and classified in the categories "annotated", "novel_acceptor", "novel_donor", "novel_combo", "novel_exon_skip", "ambig_gene" and "unannotated". Junctions were then clustered by their shared acceptor or donor sites and their counts were normalised. Coverage data associated with each junction was loaded and normalised for each sample. Junction and coverage normalised counts were scored using the z-score approach. Resulting z-scores for each patient were placed into a isolation forest model to rank splicing by how abberrant they looked in relation to controls. The rank for the most abnormal splicing event in ATG7 was then extracted. Sashimi plots describing the splicing across the ATG7 patients were plotted for each of the 5 patients.

2.1 OUTRIDER methods

2.1.1 Filter by expression

  • There are 12339 genes remaining after filtering by expression.
knitr::include_graphics(here::here("results/aberrant_expression/OUTRIDER/50_controls/FPKM.png"))

knitr::include_graphics(here::here("results/aberrant_expression/OUTRIDER/50_controls/expressed_genes.png"))

2.1.2 Pre-correction gene count heatmaps

knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_cor_pre_correction.png"))

knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_gene_sample_pre_correction.png"))

2.1.3 Optimising the autoencoder correction dimension

knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/encoding_dim.png"))

2.1.4 Post-correction gene count heatmaps

Here, I have included one example (M0920.18). In reality, there's plots like these for each ATG7 sample as I have re-run each ATG7 with the other ATG7 samples removed.

knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_cor_post_correction_M0920.18.png"))

knitr::include_graphics(here::here("results//aberrant_expression/OUTRIDER/50_controls/count_gene_sample_post_correction_M0920.18.png"))

3 Results

3.1 Obtain gene counts

3.1.1 Generate genes-only version of Ensembl v97 gtf


cd /tools
git clone https://github.com/broadinstitute/gtex-pipeline
chmod -R 775 gtex-pipeline

# generate the genes gtf
/tools/gtex-pipeline/gene_model/collapse_annotation.py \
/data/references/ensembl/gtf_gff3/v97/Homo_sapiens.GRCh38.97.gtf \
/data/references/ensembl/gtf_gff3/v97/Homo_sapiens.GRCh38.97.genes.gtf

3.1.2 RNA-SEQC

source(here::here("scripts/aberrant_expression/get_gene_count_RSE.R"))

3.2 Which pathways/genes are differentially expressed in response to disease?

3.2.1 ATG7 expression is significantly lower in family 1

We compared the expression of ATG7 measured as fragments per kilobase of exon model per million reads mapped (FPKM) against the remaining patient samples. Family 1 was observed to have lower ATG7 expression compared to the other patients. All other families looked to have similar ATG7 expression to the remaining patient samples.

plot_fpkm <- function(gene_id, gene_count_ods, ATG7_patient_metadata){
  
  
  gene_count_ods_filtered <- gene_count_ods[names(gene_count_ods) == gene_id, ]
  
  expr_to_plot <- 
    assays(gene_count_ods_filtered)[["fpkm"]] %>% 
    as.data.frame() %>% 
    gather(key = "samp_id", value = "fpkm")  %>% 
    mutate(ATG7_patient = samp_id %in% ATG7_patient_metadata[["samp_id_tidy"]]) 
  
  families <- tibble(samp_id = c("M0920.18", "M0921.18", "M1111.19", "M1716.19", "M1856.17", "S2557"), 
                     family = c(2, 2, 3, 4, 1, 1))
  
  expr_plot <- expr_to_plot %>% 
    left_join(families) %>% 
    mutate(family = ifelse(is.na(family), "non-ATG7 sample", family)) %>% 
    ggplot(aes(x = ATG7_patient, y = fpkm)) + 
    geom_boxplot() +
    geom_jitter(aes(colour = family)) + 
    scale_x_discrete(name = "ATG7 Patient") + 
    scale_y_continuous(name = "FPKM") + 
    scale_colour_manual(values = c(ggpubr::get_palette("npg", 4), "black")) +
    theme_bw()

  return(expr_plot)
  
}

ATG7_expr_plot <- plot_fpkm(gene_id = "ENSG00000197548", 
                            gene_count_ods,
                            ATG7_patient_metadata)

ATG7_expr_plot

ggsave(ATG7_expr_plot, 
       filename = "ATG7_expr_plot.png", 
       path = here::here("results/tables_figures"), 
       width = 8, height = 6, dpi = 600)

3.2.2 ATG7 is detected an expression outlier in family 1

Consistent with the expression plot above, OUTRIDER detected ATG7 as an expression outlier in family 1 with a negative z-score of -4.59.

load(here::here("results/aberrant_expression/OUTRIDER/50_controls/ATG7_res_all.rda"))

ATG7_res_ATG7_only <- ATG7_res_all %>% 
  filter(sampleID %in% ATG7_patient_metadata[["samp_id_tidy"]]) %>%
  mutate(padjust = p.adjust(pValue, method = "fdr")) %>% 
  dplyr::select(Description, everything())

ATG7_res_ATG7_only %>% 
  filter(Description == "ATG7")

3.2.3 OUTRIDER detects 9 genes that are expression outliers common to at least 2 of the ATG7 patients

OUTRIDER detected 264 unique genes as expression outliers across all samples (Table ATG7_OUTRIDER_res_signif.csv). In order to find genes that were commonly disrupted across ATG7 patients we filtered results by genes that were outliers in at least two ATG7 patient samples. The table below details the 9 genes left after this filtering. Each row represents a patient and each column represents a gene that was detected as a outlier in at least 2 of the ATG7 patients. Values in each cell represent the “z-score;pvalue” for each patient/gene respectively. Green indicates a positive z-score and red a negative z-score with respect to controls. Bold font is used for those that are significant outliers (p-value < 0.05 after FDR correction).

VPS41, HMG20A and TBC1D3L were specifically down-regulated in family 2. TBC1D3L is a GTPase activating protein for RAB5, which looks to have a role in macroautophagy (https://jcs.biologists.org/content/121/10/e1002). VPS41 has a role in vesicular transport and is suggested to be neuroprotective in AD/PD acting via an autophagy mechanism (https://www.ncbi.nlm.nih.gov/pmc/articles/PMC6276837/).

ATG7_res_ATG7_only_signif <- ATG7_res_ATG7_only %>% 
  filter(padjust <= 0.05) %>% 
  dplyr::select(Description, everything())

common_genes <- ATG7_res_ATG7_only_signif %>% 
  dplyr::count(Description) %>% 
  filter(n > 1) %>% 
  arrange(desc(n)) 

common_genes <- 
  ATG7_res_ATG7_only %>% 
  filter(Description %in% common_genes$Description, 
         sampleID %in% ATG7_patient_metadata[["samp_id_tidy"]]) %>% 
  arrange(Description) %>% 
  mutate(pval_zscore = str_c(zScore, ";", round(padjust, 2))) %>% 
  dplyr::select(`Sample ID` = sampleID, pval_zscore, gene_name = Description) %>% 
  arrange(gene_name) %>% 
  spread(key = gene_name, value = pval_zscore) %>% 
  left_join(ATG7_patient_metadata %>% 
              dplyr::select(samp_id_tidy, family), 
            by = c("Sample ID" = "samp_id_tidy")) %>% 
  dplyr::select(Family = family, everything()) %>% 
  arrange(Family)

common_genes %>% data.table::data.table()
write_delim(common_genes, here::here("results/tables_figures/OUTRIDER_common_outliers.csv"), delim = ",")

3.2.4 Exploratory analysis

3.2.4.1 Common outliers to M1856-17, M1111-19, and M1716-19

ATG7_res_ATG7_only %>% 
  filter(sampleID %in% c("M1856.17", "M1111.19", "M1716.19")) %>% 
  group_by(sampleID) %>% 
  mutate(rank = rank(padjust)) %>% 
  group_by(geneID) %>% 
  mutate(mean_rank = mean(rank)) %>% 
  ungroup() %>% 
  arrange(mean_rank) %>% 
  head(200) %>% 
  DT::datatable()
# METTL18
plot_fpkm("ENSG00000171806", gene_count_ods, ATG7_patient_metadata)

# TMEM246
plot_fpkm("ENSG00000165152", gene_count_ods, ATG7_patient_metadata)

3.2.4.2 Negative controls

There are no outlier genes detected in either of the control samples, which is reassuring.

ATG7_res_all %>% 
  filter(str_detect(sampleID, "control")) %>% 
  mutate(padjust =  p.adjust(pValue, method = "fdr")) %>% 
  filter(padjust <= 0.05)

As a double check, how many of the ATG7 outliers are common to any of the remaining patient samples.

signif_genes_remaining_samp <- ATG7_res_all %>% 
  filter(!(sampleID %in% ATG7_patient_metadata[["samp_id_tidy"]])) %>% 
  mutate(padjust =  p.adjust(pValue, method = "fdr")) %>% 
  filter(padjust <= 0.05) %>% 
  .[["Description"]]

mean(unique(ATG7_res_ATG7_only_signif$Description) %in% unique(signif_genes_remaining_samp))
## [1] 0.2083333

3.2.4.3 Gene-set enrichment

Here, we use a ranked gprofiler test with all 264 outlier genes as input. The background genes used are all genes that have been tested through OUTRIDER.

hits_gprofiler <- 
  gprofiler2::gost(query = unique(ATG7_res_ATG7_only_signif %>% 
                                    arrange(padjust) %>% 
                                    .[["geneID"]]), 
                   organism = "hsapiens", 
                   ordered_query = T, 
                   significant = T, 
                   user_threshold = 0.05, 
                   correction_method = c("g_SCS"),
                   domain_scope = c("custom"),
                   custom_bg = unique(ATG7_res_all$geneID),
                   sources = c("GO", "KEGG", "REAC"), 
                   evcodes = F)

hits_gprofiler$result %>% 
  data.table::data.table()
gprofiler2::gostplot(hits_gprofiler)

3.2.4.4 Known autophagy hits

How many known autophagy genes that are part of the KEGG pathway are found amongst the ATG7 outliers? RAB39B comes up as upregulated in one patient (M0921.18). Though this is different from the RAB9 suggested by Jack.

query_kegg_pathway <- function(pathway_name, organism = "hsa"){
  
  # get all pathways associated with humans 
  all_pathways <- keggList("pathway", organism)
  
  # filter to only selected pathway and get information about it
  selected_pathway <- all_pathways[str_detect(all_pathways, pathway_name)]
  
  print(str_c("Getting info for the pathways: ", unname(selected_pathway) %>% str_c(collapse = ", ")))
  
  
  selected_pathway_info <- keggGet(names(selected_pathway))
  
  return(selected_pathway_info)
  
}

get_pathway_gene_info <- function(kegg_pathway_info, organism = "hsa"){
  
  pathway_gene_info_all <- tibble()
  
  for(i in 1:length(kegg_pathway_info)){
    
    odds <- (1:length(kegg_pathway_info[[i]]$GENE))[(1:length(kegg_pathway_info[[i]]$GENE) %% 2) != 0]
    evens <- (1:length(kegg_pathway_info[[i]]$GENE))[(1:length(kegg_pathway_info[[i]]$GENE) %% 2) == 0]
    
    pathway_gene_info <- 
      tibble(pathway_kegg_code = kegg_pathway_info[[i]]$ENTRY, 
             pathway_name = kegg_pathway_info[[i]]$NAME, 
             organism = organism, 
             kegg_gene_entry = kegg_pathway_info[[i]]$GENE[odds], 
             gene = kegg_pathway_info[[i]]$GENE[evens]) %>% 
      separate(col = "gene", into = c("gene_name", "gene_desc"), sep = "; ")
    
    pathway_gene_info_all <- pathway_gene_info_all %>% bind_rows(pathway_gene_info)
    
  }
  
  # for(i in 1:nrow(pathway_gene_info)){ 
  #   
  #   gene_info <- keggGet(str_c(organism, ":", pathway_gene_info$kegg_gene_entry)[i])
  #   pathway_gene_info$gene_identifiers[i] <- gene_info[[1]]$DBLINKS %>% str_c(collapse = ";")
  # 
  # }

  return(pathway_gene_info_all)
  
}

autophagy_kegg_info <- query_kegg_pathway(pathway_name = "Autophagy", organism = "hsa")
## [1] "Getting info for the pathways: Autophagy - other - Homo sapiens (human), Autophagy - animal - Homo sapiens (human)"
autophagy_gene_info <- get_pathway_gene_info(kegg_pathway_info = autophagy_kegg_info, organism = "hsa")

ATG7_res_ATG7_only_signif[ATG7_res_ATG7_only_signif$Description %in% unique(autophagy_gene_info$gene_name),]
plot_fpkm(gene_id = "ENSG00000155961", gene_count_ods, ATG7_patient_metadata)

NFE2L2 may be up-regulated in ATG7 patients. NFE2L2 is not an expression outlier and just misses the p-value cut-off for DESeq2 (p-value: 0.07). However, when plotting the expression of NFE2L2 it looks like it may be up-regulated in ATG7 patients.

ATG7_res_ATG7_only %>% filter(Description == "NFE2L2")
plot_fpkm(gene_id = "ENSG00000116044", gene_count_ods, ATG7_patient_metadata)

4 Reproducibility

## ─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────
##  setting  value                       
##  version  R version 3.6.1 (2019-07-05)
##  os       Ubuntu 16.04.6 LTS          
##  system   x86_64, linux-gnu           
##  ui       X11                         
##  language (EN)                        
##  collate  en_GB.UTF-8                 
##  ctype    en_GB.UTF-8                 
##  tz       Europe/London               
##  date     2020-10-15                  
## 
## ─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────
##  package              * version  date       lib source                              
##  abind                  1.4-5    2016-07-21 [2] CRAN (R 3.6.1)                      
##  annotate               1.64.0   2019-10-29 [1] Bioconductor                        
##  AnnotationDbi        * 1.48.0   2019-10-29 [1] Bioconductor                        
##  askpass                1.1      2019-01-13 [2] CRAN (R 3.6.1)                      
##  assertthat             0.2.1    2019-03-21 [2] CRAN (R 3.6.1)                      
##  backports              1.1.10   2020-09-15 [1] CRAN (R 3.6.1)                      
##  base64enc              0.1-3    2015-07-28 [2] CRAN (R 3.6.1)                      
##  BBmisc                 1.11     2017-03-10 [1] CRAN (R 3.6.1)                      
##  Biobase              * 2.46.0   2019-10-29 [1] Bioconductor                        
##  BiocFileCache          1.10.2   2019-11-08 [1] Bioconductor                        
##  BiocGenerics         * 0.32.0   2019-10-29 [1] Bioconductor                        
##  BiocParallel         * 1.20.1   2019-12-21 [1] Bioconductor                        
##  biomaRt                2.42.1   2020-03-26 [1] Bioconductor                        
##  Biostrings             2.54.0   2019-10-29 [1] Bioconductor                        
##  bit                    4.0.4    2020-08-04 [1] CRAN (R 3.6.1)                      
##  bit64                  4.0.5    2020-08-30 [1] CRAN (R 3.6.1)                      
##  bitops                 1.0-6    2013-08-17 [2] CRAN (R 3.6.1)                      
##  blob                   1.2.1    2020-01-20 [1] CRAN (R 3.6.1)                      
##  broom                  0.5.6    2020-04-20 [1] CRAN (R 3.6.1)                      
##  car                    3.0-8    2020-05-21 [1] CRAN (R 3.6.1)                      
##  carData                3.0-4    2020-05-22 [1] CRAN (R 3.6.1)                      
##  caTools                1.18.0   2020-01-17 [1] CRAN (R 3.6.1)                      
##  cellranger             1.1.0    2016-07-27 [2] CRAN (R 3.6.1)                      
##  checkmate              2.0.0    2020-02-06 [1] CRAN (R 3.6.1)                      
##  cli                    2.1.0    2020-10-12 [1] CRAN (R 3.6.1)                      
##  cluster                2.1.0    2019-06-19 [4] CRAN (R 3.6.1)                      
##  codetools              0.2-16   2018-12-24 [4] CRAN (R 3.6.1)                      
##  colorspace             1.4-1    2019-03-18 [2] CRAN (R 3.6.1)                      
##  crayon                 1.3.4    2017-09-16 [2] CRAN (R 3.6.1)                      
##  crosstalk              1.1.0.1  2020-03-13 [1] CRAN (R 3.6.1)                      
##  curl                   4.3      2019-12-02 [1] CRAN (R 3.6.1)                      
##  data.table           * 1.13.0   2020-07-24 [1] CRAN (R 3.6.1)                      
##  DBI                    1.1.0    2019-12-15 [2] CRAN (R 3.6.1)                      
##  dbplyr                 1.4.4    2020-05-27 [1] CRAN (R 3.6.1)                      
##  DelayedArray         * 0.12.3   2020-04-09 [1] Bioconductor                        
##  dendextend             1.14.0   2020-08-26 [1] CRAN (R 3.6.1)                      
##  DESeq2                 1.26.0   2019-10-29 [1] Bioconductor                        
##  digest                 0.6.25   2020-02-23 [1] CRAN (R 3.6.1)                      
##  dplyr                * 1.0.2    2020-08-18 [1] CRAN (R 3.6.1)                      
##  DT                     0.14     2020-06-24 [1] CRAN (R 3.6.1)                      
##  ellipsis               0.3.1    2020-05-15 [1] CRAN (R 3.6.1)                      
##  evaluate               0.14     2019-05-28 [2] CRAN (R 3.6.1)                      
##  fansi                  0.4.1    2020-01-08 [1] CRAN (R 3.6.1)                      
##  farver                 2.0.3    2020-01-16 [1] CRAN (R 3.6.1)                      
##  fastmap                1.0.1    2019-10-08 [1] CRAN (R 3.6.1)                      
##  forcats              * 0.5.0    2020-03-01 [2] CRAN (R 3.6.1)                      
##  foreach                1.5.0    2020-03-30 [1] CRAN (R 3.6.1)                      
##  foreign                0.8-72   2019-08-02 [4] CRAN (R 3.6.1)                      
##  Formula                1.2-3    2018-05-03 [2] CRAN (R 3.6.1)                      
##  fs                     1.5.0    2020-07-31 [1] CRAN (R 3.6.1)                      
##  gclus                  1.3.2    2019-01-07 [1] CRAN (R 3.6.1)                      
##  gdata                  2.18.0   2017-06-06 [1] CRAN (R 3.6.1)                      
##  genefilter             1.68.0   2019-10-29 [1] Bioconductor                        
##  geneplotter            1.64.0   2019-10-29 [1] Bioconductor                        
##  generics               0.0.2    2018-11-29 [1] CRAN (R 3.6.1)                      
##  GenomeInfoDb         * 1.22.1   2020-03-27 [1] Bioconductor                        
##  GenomeInfoDbData       1.2.2    2020-08-13 [1] Bioconductor                        
##  GenomicAlignments      1.22.1   2019-11-12 [1] Bioconductor                        
##  GenomicFeatures      * 1.38.2   2020-02-15 [1] Bioconductor                        
##  GenomicRanges        * 1.38.0   2019-10-29 [1] Bioconductor                        
##  ggplot2              * 3.3.2    2020-06-19 [1] CRAN (R 3.6.1)                      
##  ggpubr                 0.4.0    2020-06-27 [1] CRAN (R 3.6.1)                      
##  ggsci                  2.9      2018-05-14 [1] CRAN (R 3.6.1)                      
##  ggsignif               0.6.0    2019-08-08 [1] CRAN (R 3.6.1)                      
##  glue                   1.4.2    2020-08-27 [1] CRAN (R 3.6.1)                      
##  gplots                 3.0.4    2020-07-05 [1] CRAN (R 3.6.1)                      
##  gprofiler2             0.1.9    2020-04-23 [1] CRAN (R 3.6.1)                      
##  gridExtra              2.3      2017-09-09 [2] CRAN (R 3.6.1)                      
##  gtable                 0.3.0    2019-03-25 [2] CRAN (R 3.6.1)                      
##  gtools                 3.8.2    2020-03-31 [1] CRAN (R 3.6.1)                      
##  haven                  2.3.1    2020-06-01 [1] CRAN (R 3.6.1)                      
##  heatmaply              1.1.1    2020-08-25 [1] CRAN (R 3.6.1)                      
##  here                   0.1      2017-05-28 [1] CRAN (R 3.6.1)                      
##  Hmisc                  4.4-1    2020-08-10 [1] CRAN (R 3.6.1)                      
##  hms                    0.5.3    2020-01-08 [1] CRAN (R 3.6.1)                      
##  htmlTable              2.1.0    2020-09-16 [1] CRAN (R 3.6.1)                      
##  htmltools              0.5.0    2020-06-16 [1] CRAN (R 3.6.1)                      
##  htmlwidgets            1.5.1    2019-10-08 [1] CRAN (R 3.6.1)                      
##  httpuv                 1.5.2    2019-09-11 [2] CRAN (R 3.6.1)                      
##  httr                   1.4.2    2020-07-20 [1] CRAN (R 3.6.1)                      
##  IRanges              * 2.20.2   2020-01-13 [1] Bioconductor                        
##  iterators              1.0.12   2019-07-26 [2] CRAN (R 3.6.1)                      
##  jpeg                   0.1-8.1  2019-10-24 [1] CRAN (R 3.6.1)                      
##  jsonlite               1.7.1    2020-09-07 [1] CRAN (R 3.6.1)                      
##  KEGGREST             * 1.24.1   2019-10-08 [1] Bioconductor                        
##  KernSmooth             2.23-17  2020-04-26 [1] CRAN (R 3.6.1)                      
##  knitr                  1.29     2020-06-23 [1] CRAN (R 3.6.1)                      
##  labeling               0.3      2014-08-23 [2] CRAN (R 3.6.1)                      
##  later                  1.1.0.1  2020-06-05 [1] CRAN (R 3.6.1)                      
##  lattice                0.20-38  2018-11-04 [4] CRAN (R 3.6.1)                      
##  latticeExtra           0.6-29   2019-12-19 [2] CRAN (R 3.6.1)                      
##  lazyeval               0.2.2    2019-03-15 [2] CRAN (R 3.6.1)                      
##  lifecycle              0.2.0    2020-03-06 [1] CRAN (R 3.6.1)                      
##  locfit                 1.5-9.4  2020-03-25 [1] CRAN (R 3.6.1)                      
##  lubridate              1.7.9    2020-06-08 [1] CRAN (R 3.6.1)                      
##  magrittr               1.5      2014-11-22 [2] CRAN (R 3.6.1)                      
##  MASS                   7.3-51.6 2020-04-26 [1] CRAN (R 3.6.1)                      
##  Matrix                 1.2-18   2019-11-27 [1] CRAN (R 3.6.1)                      
##  matrixStats          * 0.57.0   2020-09-25 [1] CRAN (R 3.6.1)                      
##  memoise                1.1.0    2017-04-21 [2] CRAN (R 3.6.1)                      
##  mime                   0.9      2020-02-04 [1] CRAN (R 3.6.1)                      
##  modelr                 0.1.5    2019-08-08 [2] CRAN (R 3.6.1)                      
##  munsell                0.5.0    2018-06-12 [2] CRAN (R 3.6.1)                      
##  nlme                   3.1-148  2020-05-24 [1] CRAN (R 3.6.1)                      
##  nnet                   7.3-12   2016-02-02 [4] CRAN (R 3.6.1)                      
##  openssl                1.4.3    2020-09-18 [1] CRAN (R 3.6.1)                      
##  openxlsx               4.1.5    2020-05-06 [1] CRAN (R 3.6.1)                      
##  OUTRIDER             * 1.5.1    2020-09-17 [1] Github (gagneurlab/OUTRIDER@7fbdaaf)
##  pcaMethods             1.78.0   2019-10-29 [1] Bioconductor                        
##  pheatmap               1.0.12   2019-01-04 [1] CRAN (R 3.6.1)                      
##  pillar                 1.4.6    2020-07-10 [1] CRAN (R 3.6.1)                      
##  pkgconfig              2.0.3    2019-09-22 [2] CRAN (R 3.6.1)                      
##  plotly                 4.9.2.1  2020-04-04 [1] CRAN (R 3.6.1)                      
##  plyr                   1.8.6    2020-03-03 [2] CRAN (R 3.6.1)                      
##  png                    0.1-7    2013-12-03 [1] CRAN (R 3.6.1)                      
##  prettyunits            1.1.1    2020-01-24 [1] CRAN (R 3.6.1)                      
##  progress               1.2.2    2019-05-16 [2] CRAN (R 3.6.1)                      
##  promises               1.1.1    2020-06-09 [1] CRAN (R 3.6.1)                      
##  PRROC                  1.3.1    2018-06-19 [1] CRAN (R 3.6.1)                      
##  purrr                * 0.3.4    2020-04-17 [1] CRAN (R 3.6.1)                      
##  R6                     2.4.1    2019-11-12 [1] CRAN (R 3.6.1)                      
##  rappdirs               0.3.1    2016-03-28 [1] CRAN (R 3.6.1)                      
##  RColorBrewer           1.1-2    2014-12-07 [2] CRAN (R 3.6.1)                      
##  Rcpp                   1.0.5    2020-07-06 [1] CRAN (R 3.6.1)                      
##  RCurl                  1.98-1.2 2020-04-18 [1] CRAN (R 3.6.1)                      
##  readr                * 1.4.0    2020-10-05 [1] CRAN (R 3.6.1)                      
##  readxl                 1.3.1    2019-03-13 [2] CRAN (R 3.6.1)                      
##  registry               0.5-1    2019-03-05 [1] CRAN (R 3.6.1)                      
##  reprex                 0.3.0    2019-05-16 [2] CRAN (R 3.6.1)                      
##  reshape2               1.4.4    2020-04-09 [2] CRAN (R 3.6.1)                      
##  rio                    0.5.16   2018-11-26 [1] CRAN (R 3.6.1)                      
##  rlang                  0.4.8    2020-10-08 [1] CRAN (R 3.6.1)                      
##  rmarkdown              2.3      2020-06-18 [1] CRAN (R 3.6.1)                      
##  rpart                  4.1-15   2019-04-12 [4] CRAN (R 3.6.1)                      
##  rprojroot              1.3-2    2018-01-03 [1] CRAN (R 3.6.1)                      
##  Rsamtools              2.2.3    2020-02-23 [1] Bioconductor                        
##  RSQLite                2.2.1    2020-09-30 [1] CRAN (R 3.6.1)                      
##  rstatix                0.6.0    2020-06-18 [1] CRAN (R 3.6.1)                      
##  rstudioapi             0.11     2020-02-07 [2] CRAN (R 3.6.1)                      
##  rtracklayer            1.46.0   2019-10-29 [1] Bioconductor                        
##  rvest                  0.3.5    2019-11-08 [1] CRAN (R 3.6.1)                      
##  S4Vectors            * 0.24.4   2020-04-09 [1] Bioconductor                        
##  scales                 1.1.1    2020-05-11 [1] CRAN (R 3.6.1)                      
##  seriation              1.2-8    2019-08-27 [1] CRAN (R 3.6.1)                      
##  sessioninfo          * 1.1.1    2018-11-05 [2] CRAN (R 3.6.1)                      
##  shiny                  1.5.0    2020-06-23 [1] CRAN (R 3.6.1)                      
##  stringi                1.5.3    2020-09-09 [1] CRAN (R 3.6.1)                      
##  stringr              * 1.4.0    2019-02-10 [2] CRAN (R 3.6.1)                      
##  SummarizedExperiment * 1.16.1   2019-12-19 [1] Bioconductor                        
##  survival               3.2-3    2020-06-13 [1] CRAN (R 3.6.1)                      
##  tibble               * 3.0.4    2020-10-12 [1] CRAN (R 3.6.1)                      
##  tidyr                * 1.1.2    2020-08-27 [1] CRAN (R 3.6.1)                      
##  tidyselect             1.1.0    2020-05-11 [1] CRAN (R 3.6.1)                      
##  tidyverse            * 1.3.0    2019-11-21 [1] CRAN (R 3.6.1)                      
##  TSP                    1.1-10   2020-04-17 [1] CRAN (R 3.6.1)                      
##  vctrs                  0.3.4    2020-08-29 [1] CRAN (R 3.6.1)                      
##  viridis                0.5.1    2018-03-29 [2] CRAN (R 3.6.1)                      
##  viridisLite            0.3.0    2018-02-01 [2] CRAN (R 3.6.1)                      
##  webshot                0.5.2    2019-11-22 [1] CRAN (R 3.6.1)                      
##  withr                  2.3.0    2020-09-22 [1] CRAN (R 3.6.1)                      
##  xfun                   0.18     2020-09-29 [1] CRAN (R 3.6.1)                      
##  XML                    3.99-0.3 2020-01-20 [2] CRAN (R 3.6.1)                      
##  xml2                   1.3.2    2020-04-23 [1] CRAN (R 3.6.1)                      
##  xtable                 1.8-4    2019-04-21 [1] CRAN (R 3.6.1)                      
##  XVector                0.26.0   2019-10-29 [1] Bioconductor                        
##  yaml                   2.2.1    2020-02-01 [1] CRAN (R 3.6.1)                      
##  zip                    2.0.4    2019-09-01 [1] CRAN (R 3.6.1)                      
##  zlibbioc               1.32.0   2019-10-29 [1] Bioconductor                        
## 
## [1] /home/dzhang/R/x86_64-pc-linux-gnu-library/3.6
## [2] /usr/local/lib/R/site-library
## [3] /usr/lib/R/site-library
## [4] /usr/lib/R/library